function ids = get_clids(ptstruct,whichgroup,dostats,doanalys)
% Get patient ids for class defined in the group requested
% It also performs some basic analysis

N = length(ptstruct);

switch whichgroup
    case 'mort28d'
        ids{1} = []; % Survivors
        ids{2} = []; % Non-survivors
        
        N1 = 0; % Patients with only one ICU
        N2 = 0; % Patients with only one ICU and with data
        
        for pt = 1 : N
            % Process patients with only one ICU admission
            if ptstruct(pt).nicus(1) == 1
                N1 = N1 + 1;
                % Select patients who have all data variables
                if ~isempty(ptstruct(pt).data{1}.hr) && ...
                        ~isempty(ptstruct(pt).data{1}.pulsep) && ...
                        ~isempty(ptstruct(pt).data{1}.uo) && ...
                        ~isempty(ptstruct(pt).data{1}.lactate)
                    N2 = N2 + 1;
                    nmort = ptstruct(pt).clinics.mort28d; % mortality category
                    ids{nmort+1} = [ids{nmort+1} pt];
                end
            end
            % Process patients with more than one ICU admission
            % ?
        end
        
    otherwise
        ids = [];
end

if dostats
    switch whichgroup
        case 'mort28d'
            for i = 1 : 2
                los{i} = []; age{i} = [];
                ismale{i} = 0; isvaso{i} = 0;
                for pt = ids{i} % Survivors
                    los{i} = [los{i} ptstruct(pt).los];
                    age{i} = [age{i} ptstruct(pt).age];
                    ismale{i} = ismale{i} + ptstruct(pt).ismale;
                    isvaso{i} = isvaso{i} + ptstruct(pt).clinics.vp;
                end
            end
            fprintf('\nTotal number of patients = %d\n',N);
            fprintf('Number of patients with 1 ICU adm = %d\n',N1);
            fprintf('Number of patients with 1 ICU adm and data = %d\n',N2);
            fprintf('\n Survivors (N = %d, %.1f%%):\n',length(ids{1}), ...
                length(ids{1})/(length(ids{1}) + length(ids{2}))*100);
            fprintf('   LoS: %.1f (%.1f) [%.2f - %.2f]',mean(los{1}), ...
                std(los{1}),min(los{1}),max(los{1}));
            fprintf('\n   Age: %.1f (%.1f) [%.2f - %.2f]',mean(age{1}), ...
                std(age{1}),min(age{1}),max(age{1}));
            fprintf('\n   Gender (males): %d (%.1f%%)',ismale{1}, ...
                ismale{1}/length(ids{1})*100);
            fprintf('\n   Vaso therapy: %d (%.1f%%)',isvaso{1}, ...
                isvaso{1}/length(ids{1})*100);
            fprintf('\n');
            fprintf('\n\n Non-Survivors (N = %d, %.1f%%):\n',length(ids{2}), ...
                length(ids{2})/(length(ids{1}) + length(ids{2}))*100);
            fprintf('   LoS: %.1f (%.1f) [%.2f - %.2f]',mean(los{2}), ...
                std(los{2}),min(los{2}),max(los{2}));
            fprintf('\n   Age: %.1f (%.1f) [%.2f - %.2f]',mean(age{2}), ...
                std(age{2}),min(age{2}),max(age{2}));
            fprintf('\n   Gender (males): %d (%.1f%%)',ismale{2}, ...
                ismale{2}/length(ids{2})*100);
            fprintf('\n   Vaso therapy: %d (%.1f%%)',isvaso{2}, ...
                isvaso{2}/length(ids{2})*100);
            fprintf('\n\n');
            
        otherwise
            ids = [];
    end
end

if doanalys
    dischdif{5} = []; % in hours
    for pt = ids{2}
        % Run analysis for all variables
        dischdif{1} = [dischdif{1}, ...
            ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.hr(:,1))))/60];
        dischdif{2} = [dischdif{2}, ...
            ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.pulsep(:,1))))/60];
        dischdif{3} = [dischdif{3}, ...
            ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.lactate(:,1))))/60];
        dischdif{4} = [dischdif{4}, ...
            ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.uo(:,1))))/60];
        
         dischdif{5} = [dischdif{5}, ...
             min([((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.hr(:,1))))/60, ...
             ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.pulsep(:,1))))/60, ...
             ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.lactate(:,1))))/60, ...
             ((ptstruct(pt).los*24*60) - (max(ptstruct(pt).data{1}.uo(:,1))))/60])];
    end
    figure;
    subplot(221); hist(dischdif{1},-5:1000); title('HR'); xlabel('Hours');
    subplot(222); hist(dischdif{2},-5:1000); title('PP'); xlabel('Hours');
    subplot(223); hist(dischdif{3},-5:1000); title('Lactate'); xlabel('Hours');
    subplot(224); hist(dischdif{4},-5:1000); title('UO'); xlabel('Hours');
    figure;
    hist(dischdif{5},-1:200); title('ALL'); xlabel('Hours');
end

end

